Protein-Level Heatmaps

Load packages

suppressPackageStartupMessages({
  library(tidyverse)
  library(pheatmap)
  library(here)
  library(ggiraph)
  library(ggrepel)
})

Differential Expression (DTE & DGE) Heatmaps

Load Transcript CPM & DGE data

Load gene CPM & DGE data

Define Subsets

top_tx_df <- tx_df %>%
  mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC_DEG)) %>%
  arrange(desc(is_sig), desc(abs_logFC))

top100_tx <- top_tx_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx <- head(top100_tx, 50)

top_gene_df <- gene_df %>%
  mutate(is_sig = FDR_DEG < 0.05, abs_logFC = abs(logFC)) %>%
  arrange(desc(is_sig), desc(abs_logFC))

top100_genes <- top_gene_df %>% slice_head(n = 100) %>% pull(row_id)
top50_genes <- head(top100_genes, 50)

Heatmap function for transcript CPM

plot_tx_heatmap <- function(tx_ids, title) {
  mat <- cpm_matrix[rownames(cpm_matrix) %in% tx_ids, , drop = FALSE]
  mat <- log2(mat + 1)
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  # Use the modified label map with PB accessions and asterisks
  rownames(mat) <- tx_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Heatmap function for gene CPM

plot_gene_heatmap <- function(gene_ids, title) {
  mat <- gene_cpm[rownames(gene_cpm) %in% gene_ids, , drop = FALSE]
  mat <- log2(mat + 1)
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  rownames(mat) <- gene_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Plot Top Transcript & Gene Heatmaps

plot_tx_heatmap(top50_tx, "Top 50 Differentially Expressed Transcripts")

plot_tx_heatmap(top100_tx, "Top 100 Differentially Expressed Transcripts")

plot_gene_heatmap(top50_genes, "Top 50 Differentially Expressed Genes")

plot_gene_heatmap(top100_genes, "Top 100 Differentially Expressed Genes")

Differential Transcript Usage (DTU) Heatmaps

Load DTU Data

Define DTU Subsets

top_tx_dtu_df <- tx_dtu_df %>%
  arrange(adj.p.value_DTU, desc(abs(lf_DTU)))

top100_tx_dtu <- top_tx_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_tx_dtu <- head(top100_tx_dtu, 50)

top_gene_dtu_df <- gene_dtu_df %>%
  arrange(adj_p.value_DTU, desc(abs(lf_DTU)))

top100_gene_dtu <- top_gene_dtu_df %>% slice_head(n = 100) %>% pull(row_id)
top50_gene_dtu <- head(top100_gene_dtu, 50)

Heatmap function for DTU – note, it is not common to show DTU with heatmaps, so I am throwing in a volcano plot as well

plot_dtu_heatmap <- function(mat, row_label_map, row_ids, title) {
  mat <- mat[rownames(mat) %in% row_ids, , drop = FALSE]
  mat <- mat[rowSums(mat) > 0, , drop = FALSE]
  mat <- mat[apply(mat, 1, function(x) sd(x, na.rm = TRUE) > 0), , drop = FALSE]
  if (nrow(mat) == 0) return(NULL)

  rownames(mat) <- row_label_map[rownames(mat)]

  annotation_df <- data.frame(
    sample_name = colnames(mat),
    group = ifelse(grepl("WT", colnames(mat)), "WT", "Q157R")
  ) %>% column_to_rownames("sample_name")

  annotation_colors <- list(group = c(Q157R = "#336B87", WT = "#E99787"))

  col_order <- c(grep("Q157R", colnames(mat), value = TRUE),
                 grep("WT", colnames(mat), value = TRUE))
  mat <- mat[, col_order]
  annotation_df <- annotation_df[col_order, , drop = FALSE]

  n_rows <- nrow(mat)
  plot_height <- min(10, 5 + n_rows * 0.15)
  fontsize_row <- if (n_rows <= 30) 10 else if (n_rows <= 60) 8 else 6

  pheatmap(mat, scale = "row", cluster_rows = TRUE, cluster_cols = FALSE,
           annotation_col = annotation_df, annotation_colors = annotation_colors,
           show_colnames = TRUE, show_rownames = TRUE, main = title,
           fontsize_row = fontsize_row, fontsize_col = 10, height = plot_height)
}

Plot DTU Transcript Heatmaps

plot_dtu_heatmap(frac_matrix, frac_label_map, top100_tx_dtu, "Top 100 DTU Transcripts")

plot_dtu_heatmap(frac_matrix, frac_label_map, top50_tx_dtu, "Top 50 DTU Transcripts")

Plot DTU Gene Heatmaps

plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top100_gene_dtu, "Top 100 DTU Genes")

plot_dtu_heatmap(gene_frac_matrix, gene_frac_label_map, top50_gene_dtu, "Top 50 DTU Genes")

Interactive Volcano plot